ゼロから始めるFFT 〜実装編〜

（実装編以外は）ないです。

FFT is 何

Fast Fourier Transformの略です。

Fast Fourier Transformと言っても$ \int_{\operatorname{hoge}}f(x)\mathrm{e}^{-2\pi\mathrm{i}xt}\mathrm{d}xみたいな積分計算をササッとするわけではなく、離散フーリエ変換を計算機上で行うアルゴリズムです。

離散フーリエ変換(DFT) is 何

長さ$ Nの列$ (a_i)_{i=0,1,\ldots,N-1}と$ 1の原始$ N乗根$ \zetaについて、$ (a_i)を離散フーリエ変換したものは長さ$ Nの列

$ \left(\hat{a}_j\right)_{j=0,1,\ldots,N-1}=\left(\sum_{i=0}^{N-1}a_i\zeta^{ij}\right)_{j=0,1,\ldots,N-1}

である。

シンプルだな！ガハハ

以下、$ N=2^k\ (k\text{は正整数})とします。畳み込みとかをする範囲においてはゼロ埋めでなんとかなるので、なる

なぜFFTがあるのか

$ N=1048576くらいで畳み込みがしたいことがあるかもしれないので、定義にしたがって計算を行うとだいたい$ N^2=1099511627776回くらいの演算が必要になるので、困る

こういうときに$ N\log N=20971520回くらいの演算で回ると、嬉しい

だから、FFTを開発する必要が、あったんですね

以下、Cooley–Tukey FFT algorithmについて話します。他にもいろいろなFFTのアルゴリズムがあるんじゃ(0埋めをしなくてもいいやつもあるぞい)

FFTのアイデア

$ a_{N/2}がどう足されるかを考えると、

$ \hat{a}_0=a_0+\cdots+a_{N/2}+\cdots

$ \hat{a}_1=a_0+\cdots+\zeta^{N/2}a_{N/2}+\cdots=a_0+\cdots-a_{N/2}+\cdots

$ \hat{a}_2=a_0+\cdots+\zeta^{N}a_{N/2}+\cdots=a_0+\cdots+a_{N/2}+\cdots

$ \vdots

みたいになっており、実は$ a_{N/2}は$ \pm a_{N/2}のどちらかの形でしか足されていないことがわかります。

同様に、$ a_kと$ a_{k+N/2}のペアは$ \zeta^{kj}\left(a_k\pm a_{k+N/2}\right)の形にまとめることができます（自分で確かめてみよう！）。

まとめて$ b_{k}=a_k+a_{k+N/2},b_{k+N/2}=a_k-a_{k+N/2}のようにすると、

$ \hat{a}_{2j}=\sum_{i=0}^{N/2}\zeta^{2ij}b_i

$ \hat{a}_{2j+1}=\sum_{i=0}^{N/2}\zeta^{(2j+1)i}b_{i+N/2}=\sum_{i=0}^{N/2}\zeta^{2ij}\zeta^ib_{i+N/2}

となることが確かめられます。これはよく見るとDFTの形をしている…！（観察力の鬼）

答えの偶数番目は$ (a_i+a_{i+N/2})_iのDFTであり、奇数番目は$ (\zeta^i(a_i-a_{i+N/2}))_iのDFTになっています。

ところで上で定めた$ (b_i)_iは線形時間で計算が可能でした（$ N回足したり引いたりするだけ）。

半分になりながら毎回線形時間かけるので、これは$ O(N\log N)時間です。

ご清聴ありがとうございました。

は？実装しろ

はい

上の解説を念頭に置いて次の列を処理することを考えます

a b c d e f g h

まず$ a_kと$ a_{k+N/2}を近くに寄せたいため、寄せる（この操作をbit reverseと言います理由はindexを2進表記すればわかる）

a e c g b f d h

出てきてほしい値と見比べる(z^4=-1です)

code:a

0 4 2 6 1 5 3 7

--+------------------------------------------------------

a | a + e + c + g + b + f + d + h

e | a - e + cz^2 - gz^2 + bz - fz + dz^3 - hz^3

c | a + e - c - g + bz^2 + fz^2 - dz^2 - hz^2

g | a - e - cz^2 + gz^2 + bz^3 - fz^3 + dz - hz

b | a + e + c + g - b - f - d - h

f | a - e + cz^2 - gz^2 - bz + fz - dz^3 + hz^3

d | a + e - c - g - bz^2 - fz^2 + dz^2 + hz^2

h | a - e - cz^2 + gz^2 - bz^3 + fz^3 - dz + hz

すると、結局こんな操作をすることになります。

code:a

a | a + e | a + e + c + g | ...

e | a - e | a - e + cz^2 - gz^2 | ...

c | c + g | a + e - c - g | ...

g | c - g | a - e - cz^2 + gz^2 | ...

b | b + f | b + f + d + h | ...

f | b - f | b - f + dz^2 - hz^2 | ...

d | d + h | b + f - d - h | ...

h | d - h | b - f - cz^2 + hz^2 | ...

規則はこの表とソースコードを見てエスパーしてくださいステップ$ iでは$ 2^iだけ離れたところをマージしていく感じ

code:cpp

template<typename T>

vector<T> fft(const vector<T>& a, const T& zeta){

size_t N{size(a)};

vector<T> ret(a);

for(size_t i{0}, j{0}; i < N; ++i){

if(i < j)swap(reti, retj);

size_t m{N};

while(m > 1 && j >= (m >>= 1))j -= m;

j += m;

}

size_t k{0}, tmpn{N};

while(tmpn >>= 1)++k;

stack<T> zeta_pow;

{

T zeta_tmp{zeta};

for(size_t i{0}; i < k; ++i){

zeta_pow.push(zeta_tmp);

zeta_tmp *= zeta_tmp;

}

for(size_t i{0}, j{1}; i < k; ++i, j <<= 1){

for(size_t k{0}; k < N; k += 2 * j){

T zeta_tmp{1};

for(size_t l{0}; l < j; ++l){

auto tmp = retk + l;

retk + l = tmp + retk + l + j * zeta_tmp;

retk + l + j = tmp - retk + l + j * zeta_tmp;

zeta_tmp *= zeta_pow.top();

}

zeta_pow.pop();

}

return ret;

}

#include<vector>とかusing namespace std;みたいなのは脳内で補完しておいてください

code:cpp

const auto pi{3.14159265358979323846264L};

const auto u{fft(fft(vector<complex<long double>>{1, 2, 3, 4, 5, 6, 7, 8}, std::polar<long double>(1.0, pi / 4.0)), std::polar<long double>(1.0, -pi / 4.0))};

for(const auto& i : u)cout << i / 8.0L << " ";

cout << endl;

code:output

(1,-5.42101e-20) (2,-8.13152e-20) (3,-5.42101e-20) (4,8.13152e-20) (5,5.42101e-20) (6,2.71051e-20) (7,5.42101e-20) (8,-2.71051e-20)

フーリエ変換して逆変換すると元の列に戻っているやったぜ（数値誤差）

ご清聴ありがとうございました